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Abstract 

A  simple  closure  for  filtered  reaction  of  a  reaction  progress  variable  is  analysed  in  this  study  using 
explicitly  filtered  DNS  data  of  turbulent  MILD  combustion  of  methane  for  Large  Eddy  Simulation 
(LES).  The  conditional  averages  of  major  and  minor  species  mass  fractions,  and  reaction  rate  constructed 
from  the  DNS  data  along  with  those  obtained  using  flamelet  and  Perfectly  Stirred  Reactor  (PSR)  models 
suggest  that  the  PSR  can  serve  as  a  good  canonical  reactor  for  MILD  combustion  modelling.  The  flamelet 
predictions  of  reaction  rate  are  observed  to  be  poor  because  it  does  not  include  effects  of  flame  interactions, 
which  are  abundant  in  the  MILD  combustion.  The  PSR  solution  obtained  over  a  wide  range  of  residence 
time  along  with  presumed  beta  sub-grid  PDF  seems  a  reasonable  closure  for  the  filtered  reaction  rate  for 
the  LES  filter  size  greater  than  three  flame  thermal  thicknesses.  Both  spatial  variations  and  joint  PDF  of 
modelled  and  DNS  values  of  filtered  reaction  rates  are  analysed. 

©  2014  The  Authors.  Published  by  Elsevier  Inc.  on  behalf  of  The  Combustion  Institute.  This  is  an  open 
access  article  under  the  CC  BY  license  (http://creativecommons.Org/licenses/by/3.0/). 
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1.  Introduction 

Due  to  the  need  for  combustion  systems  with 
improved  efficiency  and  reduced  emissions,  alter¬ 
native  combustion  technologies  are  being  con¬ 
stantly  explored  to  meet  these  requirements. 
Although  fuel-lean  premixed  combustion  can 
meet  these  requirements  it  is  highly  susceptible 
to  thermo-acoustic  instability.  Moderate  or 
Intense  Low-oxygen  Dilution  (MILD)  combus¬ 
tion  is  another  technology  employing  diluted 
and  preheated  reactants  which  can  overcome  the 
shortcoming  of  lean  premixed  combustion.  The 
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reactant  temperature,  Tr,  is  higher  than  reactant 
mixture’s  autoignition  temperature  Tign  and  the 
temperature  rise,  A T  =  Tp  —  Tr,  is  smaller  than 
T ig„  for  combustion  under  typical  MILD  condi¬ 
tions.  The  MILD  combustion  conditions  can  be 
achieved  using  Exhaust  or  Flue  Gas  Recirculation 
(EGR  or  FRG)  or  staged  fuel  injection  methods. 

Direct  photographs  of  MILD  combustion  from 
previous  experimental  studies  [1-3]  suggest  spatially 
uniform  and  steady  combustion.  A  possible  theo¬ 
retical  explanation  for  this  uniform  combustion 
has  been  offered  on  the  basis  of  two-stage  evolu¬ 
tion,  where  fuel-rich  combustion  in  diluted  oxidiser 
is  followed  by  homogeneous  combustion  of  par¬ 
tially  oxidised  fuel  [1],  More  theoretical  analyses 
of  methane  oxidation  under  diluted  conditions  have 
been  conducted  by  using  a  Perfectly  Stirred  Reactor 
(PSR)  model  [4,1, 5-8]  suggesting  that  combustion 
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in  diluted  conditions  may  be  seen  as  turbulent 
combustion  in  distributed  regime.  Although  laser 
thermometry  images  in  [9,2,10]  suggested  distrib¬ 
uted  regime  combustion,  OH  Planar  Laser 
Induced  Fluorescence  (PLIF)  images  in  the  same 
studies  show  the  presence  of  thin  reaction  zones. 
It  has  been  recently  pointed  out  [11]  that  frequent 
interactions  of  thin  reaction  zones  in  MILD  com¬ 
bustion  produce  conditions  similar  to  distributed 
regime  combustion. 

Modelling  of  MILD  combustion  is  important 
from  the  application’s  point  of  view.  The  Eddy 
Dissipation  Concept  (EDC)  is  extensively  used 
in  previous  modelling  studies  of  MILD  combus¬ 
tion  using  both  Reynolds  Averaged  Navier- 
Stokes  (RANS)  [12-17]  and  Large  Eddy  Simula¬ 
tion  (LES)  [18]  approaches.  The  computational 
results  from  these  studies  compare  reasonably  to 
experimental  results,  suggesting  that  the  EDC 
model  is  satisfactory.  However,  its  fundamental 
assumption  of  a  homogeneous  reactor  within  fine 
structures  of  turbulent  kinetic  energy  dissipation 
has  not  been  verified  for  MILD  combustion.  It 
is  important  to  note  at  this  juncture  that  the  dis¬ 
tributed  regime  of  combustion  does  not  result 
from  only  turbulence-flame  interaction,  but  also 
by  interactions  of  reaction  zones  [11]  in  MILD 
combustion.  Large  Eddy  Simulations  involving  a 
presumed  sub-grid  probability  density  function 
(PDF)  with  canonical  reactors  have  also  been 
explored  [19,20],  The  filtered  reaction  rate,  65,  in 
this  approach  is  written  as 


co  = 


coL{c)PA(c)  dc, 


(1) 


using  a  sub-grid  PDF,  PA(c),  for  a  reaction  pro¬ 
gress  variable  (RPV)  c  by  presuming  the  PDF 
shape  for  given  values  of  Favre  filtered  RPV  c 
and  its  sub-grid  variance  a2c.  The  reaction  rate 
coL(c)  may  be  assumed  to  depend  only  on  c  and 
is  taken  from  an  appropriate  canonical  model, 
either  a  PSR  [19]  or  a  steady  flamelet  [20],  Thus, 
precomputed  and  tabulated  values  for  65  can  be 
used  during  LES  using  look-up  table  approach. 
However,  choice  for  this  model  reactor  is  unclear 
in  the  light  of  experiments  cited  above,  which 
show  both  distributed  and  flamelets  combustion. 
Thus,  a  careful  investigation  is  warranted  to 
improve  our  understanding  of  turbulent  MILD 
combustion  and  to  guide  filtered  reaction  rate 
modelling. 

This  study  attempts  to  provide  a  guidance  for 
SGS  modelling  of  reaction  rate  in  MILD  combus¬ 
tion  using  Direct  Numerical  Simulation  (DNS). 
The  aim  is  to  develop  and  validate  a  representa¬ 
tive  canonical  reactor  to  model  65  using  Eq.  (1) 
for  turbulent  MILD  combustion.  The  organisa¬ 
tion  of  this  paper  is  as  follows.  The  DNS  method 
and  MILD  combustion  condition  investigated  in 
this  study  are  briefly  described  in  Section  2.  The 


details  of  analysis  are  presented  in  Section  3  and 
results  of  this  analysis  are  discussed  in  Section  4. 
A  priori  assessment  of  filtered  reaction  rate  model 
is  discussed  in  Section  5  by  directly  comparing 
modelled  65  with  those  extracted  from  the  DNS 
results.  The  conclusions  are  summarised  in 
Section  6. 


2.  DNS  of  MILD  combustion 

Turbulent  MILD  combustion  of  methane-air 
mixture  at  atmospheric  pressure  inside  a  cubic 
domain  has  been  simulated  [11,21]  directly  using 
a  DNS  code,  SENGA2  [22],  This  code  solves  fully 
compressible  conservation  equations  for  mass, 
momentum,  internal  energy  and  species  mass  frac¬ 
tions,  Yn  for  turbulent  reacting  flows.  The  govern¬ 
ing  equations  are  discretised  on  a  uniform  mesh 
employing  a  tenth  order  central  difference  scheme 
which  gradually  reduces  to  a  fourth  order  scheme 
near  boundaries.  The  integration  in  time  is 
achieved  using  a  third  order  Runge-Kutta 
scheme.  The  methane-air  combustion  kinetics  is 
modelled  using  a  skeletal  mechanism  including 
16  species  and  25  elementary  reactions  [23],  The 
transport  and  thermo-chemical  properties  of  the 
mixture  and  individual  species  are  temperature 
dependent  and  these  species  have  non-unity  Lewis 
numbers  as  they  have  different  mass  diffusivities. 

The  fuel  and  air  are  usually  fed  through  a  high 
momentum  jet  into  a  stream  of  exhaust  gases  to 
achieve  conditions  for  MILD  combustion.  The 
high  temperature  for  either  fuel  or  air  is  usually 
achieved  by  recovering  heat  in  the  exhaust  stream. 
The  time  available  for  air,  fuel  and  exhaust  gases 
to  mix  is  usually  short  in  a  typical  MILD  combus¬ 
tion  arrangement.  This  results  in  inhomogeneous 
mixtures  made  of  unburnt  and  burnt  gases.  These 
physical  processes  are  included  in  the  DNS 
through  a  two-stage  strategy  detailed  in  Refs. 
[11,21,24].  In  the  first  stage,  mixing  between  fresh 
and  exhaust  gases  is  simulated  directly  inside  a  peri¬ 
odic  domain  without  chemical  reactions,  and  the 
duration  of  this  DNS  is  about  one  large  eddy  turn¬ 
over  time  which  is  much  shorter  than  the  autoigni¬ 
tion  delay  time  for  the  given  mixture  conditions, 
but  long  enough  to  have  partially  premixed  fresh 
and  exhaust  gases  [24],  The  turbulence  Reynolds 
number,  Re,  =  u'l0/v  =  96.1,  where  u',l0  and  v 
are  respectively  the  turbulence  intensity,  integral 
length  scale  and  kinematic  viscosity  of  reactant 
mixture,  is  similar  to  those  in  previous  experimen¬ 
tal  studies  [25-27].  In  the  second  stage,  this  inho¬ 
mogeneous  mixture  is  used  as  the  initial  and 
inflowing  mixture  for  the  combustion  simulation. 
This  inflowing  mixture  has  u' =  16.4  m/s  and 
l0  =  1 .48  mm. 

The  volume  averaged  temperature,  Tm,  of  the 
inflowing  mixture  field  is  HOOK  while  Tr  = 
Tm±0.02Tm  resulting  from  mixing  of  reactants 
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and  exhaust  gases.  The  maximum,  and  aver¬ 
aged,  (Xo 2,r),  mole  fractions  of  oxygen  in  reactant 
mixture  are  respectively  0.035  and  0.025  which  are 
similar  to  those  for  a  typical  MILD  combustion 
condition  [5].  The  stoichiometric,  £SI,  and  aver¬ 
aged,  (£),  values  of  Bilger  mixture  fraction  are 
0.010  and  0.008  respectively  because  of  intense 
dilution.  The  root  mean  square  (RMS)  value  of 
the  mixture  fraction  fluctuation,  is  around  1% 
of  {£}.  The  Damkohler  number.  Da  =  (l0/SF)/ 
(u'/SL),  and  Karlovitz  number,  Ka  rj  (u' /SL)3^2 / 
( l0/dF )-1^2,  are  respectively  0.69  and  11.9,  where 
SL  is  the  burning  velocity  of  a  representative 
unstrained  laminar  flame  and  SF  is  the  Zeldovich 
thickness.  This  turbulent  MILD  combustion  cor¬ 
responds  to  Case  B  in  [21]  and  Case  B1  in  [24]. 
The  laminar  flame,  called  as  MIFE  (MILD  Flame 
Element)  in  [24],  has  a  reactant  mixture  composi¬ 
tion  with  major  species  mass  fraction  values  equal 
to  volume  averaged  values  in  the  inflowing  mix¬ 
ture  used  for  the  DNS.  Thus,  the  reactant  mixture 
in  MIFE  includes  the  recirculated  gas  mixture 
excluding  minor  species.  The  burning  velocity, 
thermal  thickness,  5th  =  (Tp  —  7Y)/|Vr|max,  and 
SF  of  MIFE  are  1.66  m/s,  1.29  mm  and  0.22  mm 
respectively. 

The  DNS  domain  is  cubic  of  size  Lx  x  Lyx 
L:  =  10  X  10  x  10  mm3  with  inflow  and  non¬ 
reflecting  outflow  boundaries  specified  using 
Navier-Stokes  characteristic  boundary  condition 
[28]  in  x  and  periodic  conditions  in  y  and  z  direc¬ 
tions.  This  computational  domain  is  discretised 
using  Nx  x  Ny  x  Nz  =  384  x  384  x  384  grid  points. 
These  grid  points  ensured  that  there  are  about  25 
grid  points  inside  the  thermal  thickness,  5,h,  and  2 
mesh  points  inside  the  Kolmogorov  length  scale. 
The  simulation  was  run  for  1.5  flow-through  time, 
xD,  before  collecting  data  for  statistical  analysis. 
The  simulation  was  then  continued  for  one 
additional  flow-through  time  and  80  data  sets 
were  collected  for  detailed  analyses.  Elaborate 
discussion  of  this  DNS  data  can  be  found  in 
Ref.  [21,24],  The  DNS  results  are  investigated 
using  instantaneous  and  filtered  fields  in  the  fol¬ 
lowing  sections.  An  investigation  regarding  RANS 
modelling  of  this  DNS  data  is  reported  in  [29], 


3.  LES  filtering  operation  and  presumed  PDF 

For  the  present  analysis,  a  quantity  /(x) 
obtained  from  the  DNS  data  is  filtered  using 

7(x)  =  //(x')G(x  -  x';  A)  dx1,  (2) 


where  x  is  a  spatial  vector  and  G  is  the  Gaussian 
kernel  with  filter  size  A  given  by 


G(x  —  x';  A)  = 


nA2 


1/2 


exp 


— 6(x  —  x')2 


(3) 


The  Favre  filtered  quantity  is  obtained  using 
/  =  p//p  and  the  variance  of  Favre  fluctuation 

is?5  =p(/-7)2/p- 

The  presumed  PDF,  where  £  is  the 

sample  space  variable  for  c,  is  computed  using 

the  /^-function  with  given  values  of  c  and  c"2  as: 


Pf{  0  = 


r’a-cr 


h 


r'o-cr1 


r(r  +  s)  ’ 


c2  ( 1  —  c)  „  r(\  —  c) 


(4) 

(5) 

(6) 


The  values  of  c  and  c"1  are  obtained  from  the 
DNS  data  using  Eq.  (2)  in  this  study,  and  these 
values  would  be  obtained  by  solving  their  trans¬ 
port  equations  in  actual  LES. 


4.  Results  and  discussion 

4.1.  DNS  results 

The  RPV  is  defined  as  cT  =  (T  —  Tr)/(TP  —  Tr) 
for  this  study  and  thus  the  reaction  rate  is  given  by 
c oCT  —  Q/cp(Tp  —  Tr),  where  Q  and  cp  are  the  heat 
release  rate  and  specific  heat  capacity  respectively. 
A  typical  spatial  variation  of  this  reaction  rate 
normalised  using  prSL/dth,  where  pr  is  the  reactant 
mixture  density,  in  the  middle  x-y  plane  is  shown 
in  Fig.  1  for  t=  1.5t0.  The  superscript  “+”  here 
and  in  the  following  discussion  denotes  norma¬ 
lised  quantities  using  laminar  flame  values  appro¬ 
priately.  For  example,  length  and  reaction  rate  are 
normalised  using  1  /dlh  and  prSL/dth.  In  Fig.  1, 


x  + 

Fig.  1 .  Typical  variation  of  normalised  reaction  rate  co+ 
in  the  middle  x-y  plane  at  t=\.5xD.  Contours  of 
wCT  =  0.3,  0.5  and  0.7  are  shown  for  clarity,  where 
mCT  =  racr/foCr,max  and  ojCT  m.ix  is  the  maximum  reaction 
rate  in  the  domain. 
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regions  having  co+  >  1  occupy  a  substantial  part 
of  the  computational  volume.  Although  intense 
reaction  zones  (see  mCT  =  a>CT/coCTimnL  =  0.5  con¬ 
tour)  seem  to  have  thin  structures  locally,  these 
reaction  zones  are  highly  convoluted  resulting  in 
their  inevitable  interactions  and  such  interactions 
yield  distributed  reaction  zones  combustion 
([11]).  The  contour  of  several  oo+  levels  show  a 
patchy  appearance  of  reaction  zones  spread 
throughout  the  domain.  Similar  behaviour  of 
reaction  zones  in  MILD  combustion  may  be  per¬ 
ceived  using  OH  Planar  Laser  Induced  Fluores¬ 
cence  (PLIF)  images  in  [10]. 

The  PDF  of  cT,  in  RANS  context,  constructed 
using  samples  collected  in  y-z  plane  at  a  given  x 
location  and  from  the  entire  sampling  period  is 
shown  in  Fig.  2  for  several  x+  =  x/5,h  locations. 
Although  this  PDF  shows  a  large  probability  of 
finding  unburnt  gases  at  x+  =  0.30,  its  variation 
is  generally  broad.  Especially,  the  PDF  at 
x+  =4.1  has  a  plateau  suggesting  a  relatively  uni¬ 
form  temperature  distribution  without  strong  sca¬ 
lar  gradients  [1-3],  The  PDF  peak  increases  and 
shifts  towards  larger  cT  as  seen  in  Fig.  2  if  one 
moves  further  downstream  through  the  computa¬ 
tional  domain.  Also,  the  PDF  remains  broad  with 
non-zero  probability  for  0.5  ^  cT  <  1.  The  bimo- 
dal  behaviour  for  this  PDF  is  not  observed  imply¬ 
ing  that  the  use  of  classical  flamelets  to  model 
mean  reaction  rate  in  MILD  combustion  may 
not  be  fully  acceptable  and  alternative  approaches 
are  to  be  found.  It  is  worth  noting  here  that  this 
non-flamelet  like  PDF  can  also  result  from  fre¬ 
quent  and  extensive  interactions  of  thin  flames, 
producing  an  intense  reaction  rate  but  a  decrease 
in  the  local  scalar  gradient  [11],  Thus,  the  influ¬ 
ence  of  interactions  on  flamelets  must  be  taken 
into  account. 

The  sub-grid  PDF,  PA(cT),  constructed  using 
samples  collected  within  the  filter  volume  of  A5 
using  the  DNS  data  at  t  =  \.5zD  is  shown  in 
Fig.  3.  The  centre  of  this  volume  is  used  to 


Fig.  2.  The  PDF  of  cT  is  shown  for  several  streamwise 
(x)  locations.  This  RANS  PDF  is  constructed  using  data 
collected  over  entire  sampling  period. 


Fig.  3.  The  sub-grid  PDF,  Pa(ct),  for  three  streamwise 
positions  x+  =  1.9,  3.9  and  5.8.  Thick  grey  lines  are 
explicitly  filtered  DNS  values  and  thin  lines  are  modelled 
values.  Solid  line:  A  =  3 6lh,  and  dashed  line:  A  =  <5,/, 
(shown  only  for  x+  =  3.9). 


indicate  location  of  this  PDF  in  this  figure  and 
three  such  locations  are  noted.  Two  filter  sizes, 
A  =  S,/,  and  3 <5,*,  are  used.  The  result  for 
x+  =  3.9  is  shown  for  both  of  these  two  filter  sizes 
and  it  is  shown  only  for  A  =  3<5,/,  for  the  other  two 
locations  in  Fig.  3.  The  thick  lines  represent  the 
DNS  results  and  the  thin  lines  are  modelled 
PDFs,  Pp(cT),  calculated  using  the  f)  function  in 

Eqs.  (4)-(6)  for  the  values  of  cT  and  cf  obtained 
from  the  DNS  data.  The  sub-grid  PDFs  also  show 
neither  bimodality  nor  uniform  distribution  and 
the  non-smooth  variation  of  PA  is  because  of 
limited  samples  from  one  snap-shot  of  the  DNS 
data  used  for  the  analysis.  This  specific  compari¬ 
son  is  shown  here  so  that  a  “true’'  behaviour  at 
a  given  time  (in  a  single  snap-shot)  can  be  studied 
without  the  influence  of  time  averaging. 

The  PDF  is  narrow  for  A  =  <5,/,  and  it  is  broad 
for  the  wider,  A  =  3  5th  filter  used.  These  varia¬ 
tions  are  represented  reasonably  by  the  modelled 
PDF,  Pp.  The  sub-grid  PDFs  obtained  using  the 
DNS  data  collected  over  the  entire  sampling 
period  are  used  to  verify  the  closure  model  in 
Eq.  (1)  for  the  filtered  reaction  rate  in  the  discus¬ 
sion  below. 

4.2.  Perfectly  Stirred  Reactor  model 

Both  the  RANS  and  sub-grid  PDFs  studied  in 
the  previous  section  suggest  that  the  MILD  com¬ 
bustion  conditions  considered  in  this  study 
involve  distributed  reaction  zones  regime  combus¬ 
tion.  This  regime  for  MILD  combustion  is  not 
because  of  idealised  homogeneous  reactive  mix¬ 
ture  but  due  to  interacting  flamelets  having  com¬ 
mon  features  of  conventional  combustion  and 
some  distinctive  characteristics  of  their  own,  for 
example  the  interplay  between  autoignition  and 
flame  propagation  [21],  These  physical  processes 
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result  in  small  scalar  gradients  which  are  not  the 
usual  signature  of  typical  flamelet  combustion 
[21].  Based  on  these  observations,  a  PSR  having 
a  residence  time  representative  of  local  fluid 
dynamic  conditions  is  explored  here.  Although  a 
similar  model  was  used  in  a  previous  LES  study 
[19],  its  direct  assessment  has  not  been  conducted 
hitherto  and  thus  it  is  investigated  in  this  study. 
The  initial  diluted  reactant  mixture  composition 
is  based  on  volume  averaged  mass  fraction  values 
for  all  species  present  in  the  initial/inflowing  mix¬ 
ture  field  used  for  the  DNS.  Note  that  the  PSR 
mixture  composition  is  not  exactly  the  same  as 
that  of  MIFE  for  which  only  the  minor  species 
are  excluded  to  avoid  reactions  at  the  inflow 
boundary  [21].  Since  this  PSR  is  modelled  as  an 
unsteady  homogeneous  reactor  and  is  very 
small,  the  only  independent  variable  time,  t,  is 
converted  to  RPV  through  cT(t)  =  ( T(t )  —  T,.)/ 
{Tp  —  Tr),  where  T(t)  is  the  temperature  of  the 
reactor  at  time  t  and  Tp  is  the  reactor  temperature 
when  the  residence  time  is  very  large.  This  zero 
dimensional  unsteady  PSR  governing  equations 
with  initial  conditions  Tr  =  Tm  and  F,  =  ( Yir) 
(including  recirculated  gases  containing  radicals) 
are  integrated  in  time  using  COSILAB  [30]  until 
T(t)  reaches  a  steady  value.  Here,  the  residence 
time  is  equal  to  the  mean  convection  time  zD  from 
the  inflow  to  outflow  boundaries  in  the  DNS. 
Note  that  residence  times  ranging  from  td/4  to 
td  are  also  considered  but  it  does  not  change  the 
cT  —  coCT  relationship  unduly,  although  maximum 
value  of  cT  becomes  lower  with  a  decrease  in  res¬ 
idence  time.  Also,  it  is  found  that  the  temperature 
of  the  PSR  with  a  residence  time  of  xD  is  similar  to 
Tp  despite  a  difference  in  minor  species  composi¬ 
tion  between  MIFE  and  PSR.  The  same  chemical 
kinetics  employed  for  the  DNS  is  also  used  for  the 
PSR  calculations.  The  ignition  delay  time  defined 
based  on  the  instant  of  maximum  temperature  rise 
is  about  0.35  ms  and  this  relatively  short  delay 
time  is  due  to  the  presence  of  radicals  and  inter¬ 
mediates  in  the  initial  mixture. 

The  variations  of  F,  and  o>+  with  cT  are  shown 
in  Fig.  4.  The  scattered  data  is  from  the  DNS 
result  taken  at  t=  1.5td  from  the  entire  DNS 
domain  and  it  shows  a  larger  scatter  for  F,  and 
coCT  near  cT  ~  0.  The  range  for  this  scatter 
decreases  as  cT  increases  due  to  the  progress  of 
turbulent  mixing,  molecular  diffusion  and  chemi¬ 
cal  reactions.  The  general  trend  for  the  variation 
of  DNS  quantities  with  cT  shown  in  Fig.  4  is  rep¬ 
resented  by  the  conditional  averages,  (F,|cr)  and 
(co+ |cr)  shown  in  this  figure.  These  conditional 
averages  are  computed  using  samples  collected 
over  the  entire  sampling  period.  The  values  of 
(FcH4|cr)  decreases  as  the  temperature  increases. 
The  conditional  OH  mass  fraction  is  non-zero 
even  for  small  cT  values,  which  results  from  the 
mixing  of  OH  present  in  the  exhaust  gases  with 
reactant.  This  leads  to  non-zero  (co+| cT)  even 


(a)  (b)  (c) 


Fig.  4.  Variations  of  (a)  7Ch4  ( x  103),  (b)  F0h  ( x  103) 
and  (c)  with  eT.  The  scattered  data  is  the  DNS  result 
at  t  =  To,  thick  solid  line  is  the  conditional  averages, 
(F,|cr)  (a  and  b)  and  (co+  \ct)  (c),  obtained  from  the 
DNS  data  collected  over  the  entire  sampling  period,  red 
dashed  line  is  the  flamelet  (MIFE)  solution,  and  thin 
solid  line  is  the  PSR  result.  (For  interpretation  of  the 
references  to  color  in  this  figure  legend,  the  reader  is 
referred  to  the  web  version  of  this  article.) 

for  cT  ~  0  suggesting  that  local  mixtures  start  to 
react  as  soon  as  they  enter  the  computational 
domain.  This  is  because  of  high  Tr  and  presence 
of  chemically  active  radicals  near  the  inlet,  as 
one  would  expect  in  MILD  combustion  because 
of  recirculation  of  exhaust  gases.  These  reaction 
zones  near  the  inlet  are  spatially  non-uniform, 
{(°ZT\cr)  values  in  the  DNS  are  smaller  than  the 
PSR  value  for  cT  ~  0.  However,  there  are  some 
samples  in  the  DNS  with  as  large  co+  as  in  the 
PSR  near  cT  ~  0. 

The  results  for  the  two  models  involving  flam- 
elets  (MIFE)  and  PSR  are  compared  with  the 
DNS  results  in  Fig.  4.  Although  the  MIFE  solu¬ 
tion  represents  the  variations  of  FCH4  well  the 
comparisons  shown  for  F0h  and  reaction  rate 
are  not  good  as  in  Fig.  4.  Similar  comparisons 
are  observed  for  other  major  and  minor  species 
mass  fractions  values  (not  shown).  The  disagree¬ 
ment  seen  for  the  minor  species  and  reaction  rate 
is  because  of  the  presence  of  radicals  and  interme¬ 
diates  in  the  inflowing  mixture  field  for  the  DNS. 
These  species  cannot  be  included  for  the  MIFE, 
because  the  flame  speed  eigenvalue  does  not  exist 
when  the  radical  and  intermediates  are  present  in 
the  reactant  mixture  at  high  temperature  such  as 
in  MILD  combustion  [31],  Although  the  inflowing 
mixture  is  non-uniform,  the  zero-dimensional 
PSR  model  gives  reasonable  agreement  with  the 
conditional  averages  of  both  the  major  and  minor 
species  mass  fraction  values,  and  reaction  rate 
(only  FCh4  and  F0h  are  shown  in  Fig.  4).  The 
reaction  rate  for  low  cT  values  is  significantly  large 
for  the  PSR  model  compared  to  the  DNS  values. 
This  significant  difference  is  due  to  the  boundary 
conditions;  radicals  are  confined  to  a  relatively 
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small  portion  creating  non-uniform  reaction  rate 
field  for  low  cT  value  which  occurs  near  the  inflow 
boundary  of  the  DNS  domain.  However,  good 
comparisons  seen  in  Fig.  4  suggest  that  the  PSR 
model  is  a  suitable  candidate  for  the  MILD  condi¬ 
tions  investigated  in  this  study. 

5.  A  priori  assessment 

The  filtered  reaction  rate  of  cT  is  modelled 
using  Eq.  (1)  and  the  p  PDF  in  Eqs.  (4)-(6).  This 
modelled  value  is  denoted  as  55+ model  in  the  fol¬ 
lowing  discussion.  This  modelled  reaction  rate 
field  is  compared  to  the  filtered  reaction  rate  a>+ 
obtained  by  explicitly  filtering  the  DNS  results 

using  Eq.  (2).  The  values  of  c  and  c"2  required 
for  the  p  PDF  model  are  obtained  by  filtering 
the  DNS  values  appropriately.  Three  values  of  fil¬ 
ter  widths,  A  =  <5,/,,  2<5,/,  and  3Slh,  are  used  for  the 
assessment  discussed  in  this  subsection. 

Figure  5  shows  a  comparison  of  co+  and 
mj,  model  'n  middle  x-y  plane  at  t=  1.5td  for 


A  =  8,h,25,h  and  3S,h.  Note  that  this  plane  is  the 
same  as  in  Fig.  1,  but  only  a  small  part, 
2  ^  y+  ^  6,  is  shown  for  a  better  comparison. 
Small  convoluted  islands  of  intense  reaction  rates 
observed  in  Fig.  1  become  smoothed  due  to  filter¬ 
ing  operation  to  construct  cb+ .  The  smoothing 
becomes  more  obvious  as  the  filter  size  increases, 
which  yields  an  appearance  of  distributed  reac¬ 
tion  zones  with  co+  5=  1  spreading  over  a  large 
portion  of  the  domain.  The  modelled  reaction  rate 
co+  model  seems  to  reproduce  the  filtered  reaction 
rate  field  reasonably  well  as  shown  in  Fig.  5.  How¬ 
ever,  the  model  underestimates  the  filtered  reac¬ 
tion  rate  near  the  inflow  boundary. 

The  performance  of  this  reaction  rate  model  in 
a  broader  context  can  be  examined  by  studying 
the  joint  PDF  (JPDF)  of  cb+  and  ro+>modcl.  This 
JPDF  is  shown  in  Fig.  6  for  three  filter  widths. 
If  there  is  a  perfect  agreement  between  the  mod¬ 
elled  and  DNS  values  of  the  filtered  reaction  rate 
then  this  JPDF  will  have  a  spread  around  the 
diagonal  lines  shown  in  this  figure.  As  one  can 
observe  in  Fig.  6,  there  is  an  increased  tendency 


(a) 

Fig.  5.  Comparison  of  aj+  (top)  and  5j+ modd  (bottom)  in  the  middle  x-y  plane  at  t  =  1.5td  for  (a)  A  =  5,k,  (b)  2<5,*  and 
(c)  36,h.  Note  that  this  plane  is  the  same  as  in  Fig.  1,  but  only  2  ^  6  is  shown. 
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(a) 


(b)  (c) 


Fig.  6.  The  JPDF  of  eo+  and  co^  modd  computed  using  the  DNS  data  from  the  entire  sampling  period  for  (a)  A  =  5,k,  (b) 
25, h  and  (c)  35, h . 
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for  this  spread  for  larger  filter  widths.  The  reason 
for  this  is  as  follows.  The  conditions  of  the  current 
MILD  combustion  case  include  both  thin  and  dis¬ 
tributed  reaction  zones  as  seen  in  Fig.  1 .  Most  of 
these  thin  reaction  zones  are  filtered  out  for 
A  =  2S,i,  and  3(5,*,  while  they  are  likely  to  remain 
for  A  =  3,h.  Obviously,  the  PSR  model  does  not 
cater  for  thin  reaction  zones  resulting  in  a  poor 
comparison  for  A  =  3„,  as  seen  in  Fig.  6a.  These 
results  suggest  that  the  PSR  model  with  a  stan¬ 
dard  presumed  PDF  approach  works  reasonably 
well  for  MILD  combustion  modelling.  The  classi¬ 
cal  flamelets  model  may  have  to  be  extended  to 
include  effects  of  flame  interactions  before  it  could 
be  used  for  turbulent  MILD  combustion. 


6.  Summary 

DNS  data  of  turbulent  MILD  combustion  of 
methane  are  analysed  to  construct  a  simple  model 
for  filtered  reaction  rate  required  for  Large  Eddy 
Simulations.  This  analysis  showed  that  reaction 
zones  in  MILD  combustion  are  spread  over  large 
portion  of  the  computational  domain  and  are 
partly  thin  structures  having  some  common  fea¬ 
tures  of  conventional  premixed  combustion. 
These  thin  structures  interact  frequently  giving 
an  appearance  of  distributed  regime  combustion. 
These  views  expressed  in  an  earlier  study  [21]  are 
supported  further  by  the  marginal  PDF  of  norma¬ 
lised  temperature,  cT.  Both  sub-grid  and  RANS 
PDFs  are  analysed  and  (i  function  PDF  modelled 
using  Favre  filtered  cT  and  its  sub-grid  variance 
obtained  from  the  DNS  data  is  found  to  be  a  rea¬ 
sonable  model. 

The  conditionally  averaged  species  mass  frac¬ 
tions  and  reaction  rate  for  cT  constructed  from 
the  DNS  data  are  compared  to  those  obtained 
from  flamelet  and  PSR  models.  A  flamelet,  called 
as  MIFE,  having  the  reactant  mixture  composi¬ 
tion  corresponding  to  the  DNS  mixture  is  used. 
The  residence  time  for  the  unsteady  PSR  model 
has  a  broad  range  to  cover  all  possible  reacting 
state  that  could  exist  in  turbulent  MILD  combus¬ 
tion  for  thermochemical  and  turbulence  condi¬ 
tions  considered  here.  Although  the  flamelet 
model  compares  well  with  DNS  data  for  major 
species  mass  fractions,  the  conditional  averages 
of  minor  species  mass  fractions  and  reaction  rates 
do  not  compare  well.  This  could  be  due  to  the 
presence  of  flame  interactions  which  are  not 
included  in  the  flamelets  model.  The  avenues  to 
include  effects  of  interaction  in  flamelets  based 
approach  are  unclear  and  it  is  an  open  question 
at  this  time.  However,  it  would  also  be  interesting 
to  compare  turbulent  MILD  reaction  zones  with 
canonical  reaction  zones  other  than  MIFE  and 
PSR,  such  as  Homogeneous  Charge  Diffusion 
Ignition  (HCDI)  [32],  The  PSR  results  compare 
very  well  for  major  and  minor  species  mass 


fractions,  and  reaction  rate.  Thus,  this  model 
along  with  the  presumed  /?  sub-grid  PDF  is  used 
to  find  a  simple  closure  for  filtered  reaction  rate. 
This  modelled  reaction  rate  field  and  that 
obtained  by  filtering  the  DNS  result  agree  reason¬ 
ably  well  for  filter  size  of  A  ^  3(5,*. 
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